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Abstract 

The two-sided Rayleigh quotient iteration proposed by Ostrowski computes a pair of 
corresponding left-right eigenvectors of a matrix C . We propose a Grassmannian ver- 
sion of this iteration, i.e., its iterates are pairs of p-dimensional subspaces instead of 
one-dimensional subspaces in the classical case. The new iteration generically converges 
locally cubically to the pairs of left-right p-dimcnsional invariant subspaces of C . More- 
over, Grassmannian versions of the Rayleigh quotient iteration arc given for the general- 
ized Hermitian cigenproblcm, the Hamiltonian cigcnproblcm and the skew-Hamiltonian 
cigenproblcm. 

Keywords. Block Rayleigh quotient iteration, two-sided iteration, Grassmann manifold, 
generalized eigenproblem, Hamiltonian eigenproblem. 

AMS subject classification. 65F15 

1 Introduction 

The Rayleigh quotient iteration (RQI) is a classical method for computing eigenvectors of a 
Hermitian matrix A = A H [Par 7A\ IPar98j . The RQI is a particular inverse iteration |Ips97| 
where the shift is the Rayleigh quotient evaluated at the current iterate. The Rayleigh 
quotient is an efficient shift because in a neighborhood of any eigenvector of A it yields a 
quadratic approximation of the corresponding eigenvalue. This remarkable property endows 
the iteration with cubic rate of convergence to the eigenvectors of A (see [Par 74} IPar98j or 
the sketch of proof in |AMSV02j ). Thanks to its fast convergence, the RQI is particularly 
efficient for refining estimates of eigenvectors. 

In some cases, one has to refine an estimate of a p-dimensional invariant subspace (or 
eigenspace) of A. A reason for considering an eigenspace instead of individual eigenvec- 
tors may be that the eigenvectors themselves are ill-conditioned while the subspace is not 
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(see e.g. |Ste73j ) or just because they are not relevant for the application. Several meth- 
ods have been proposed for refining invariant subspace estimates. A quadratically conver- 
gent iteration for refining eigenspaces of arbitrary (possibly non-Hermitian) matrices was 
proposed by Chatelin [Cha844 IDem87j . It uses Newton's method for solving the Riccati 
equation obtained by expressing the eigenproblem in inhomogeneous coordinates. Similar 
Newton-like iterations for eigenspace refinement were obtained using a differential-geometric 
approach [EAS981 lLEfl2l IAMS04] : see [AMS07| for an overview. In [Smi97HSM5V02| . it was 
shown that the RQI, originally defined on the set of one-dimensional subspaces of M n , can be 
generalized to operate on the set of p-dimensional subspaces of M n . The generalized iteration, 
called block- RQI or Grassmann-RQI (because the set of the p-dimensional subspaces of M n is 
termed a Grassmann manifold) converges locally cubically to the p-dimensional eigenspaces 
of A = A H . 

It is natural to ask whether the Grassmann-RQI method can be adapted to deal with 
non-Hermitian matrices. This is the topic of the present paper. 

The underlying idea comes from Ostrowski's series of papers dedicated to the RQI |Ost59aJ . 
Let C be a nonnormal matrix. Then the quadratic approximation property of the Rayleigh 
quotient is lost (and moreover the global convergence properties of the RQI become weaker, 
see |BS90j ) . This drawback was avoided by Ostrowski jOst59bl IPar74j by considering the bi- 
lateral Rayleigh quotient p(yL, Vr) '■= y^Cynjy^yR which displays the quadratic property in 
the neighborhood of the pairs of left-right nondefective eigenvectors of C. Using this Rayleigh 
quotient as a shift, he derived a two-sided iteration (see Algorithm 12.61 below) that operates 
on pairs of vectors (or pairs of one-dimensional subspaces, since the norm is irrelevant) and 
aims at converging to pairs of left-right eigenvectors of C. The rate of convergence is in cubic 
in nondegenerate cases. The possibility of solving the two-sided RQI equations approximately 
was investigated in |HS03j . 

In the present paper, we generalize Ostrowski's two-sided RQI to operate on pairs of p- 
dimensional subspaces (instead of one-dimensional subspaces in the original iteration). The 
new iteration, called Two-Sided Grassmann-RQI (2sGRQI), converges locally cubically to 
the pairs of left-right p-dimensional eigenspaces of C (see Section [5]). Comparison between 
Chatelin's iteration and the 2sGRQI (Section [6]) shows that each method has its advantages 
and drawbacks. Main advantages of the 2sGRQI over Chatelin's iteration are the higher 
rate of convergence, the simultaneous computation of left and right eigenspaces, and the 
simpler structure of the Sylvester equations. On the other hand, the 2sGRQI does not behave 
satisfactorily when C is defective and it involves two Sylvester equations instead of one. We 
also show that in some structured eigenproblems, namely .E-(skew-)Hermitian matrices with 
E = ±E , a relation 3^z = Ey^ between left and right subspaces is invariant by the 2sGRQI 
mapping (Section [7]). In particular, this observation yields a modified one-sided Grassmann- 
RQI for the Hamiltonian eigenproblem. We report on numerical experiments in Section [8] and 
conclusions are drawn in Section [9) 

2 Preliminaries 

This paper uses a few elementary concepts related to the algebraic eigenvalue problem, such 
as principal vectors, Jordan blocks and nonlinear elementary divisors. A classical reference 
is |Wil65] . Notions of subspaces and distance between them can be found in |Ste73j . 
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The superscript H denotes the conjugate transpose. In accordance with Parlett's conven- 
tions [Par 741 IPar98j , we try to reserve the letter A for Hermitian matrices while C may denote 
any matrix. We use Grass(p, n) to denote the Grassmann manifold of the p-dimensional sub- 
spaces of C n , P n_1 to denote the projective space (i.e., the set of all one-dimensional subspaces 
of C n ), and C™ xp to denote the noncompact Stiefel manifold, i.e., the set of n-hj-p matrices 
with full rank. The space spanned by the columns of Y G C™ xp is denoted by \Y\ and called 
the span of Y. The norm of a vector x is ||a;|| = V x H x. The spectral norm of a matrix T, 
denoted by ||T||, is the largest singular value of T. The Hermitian angle Z(x,y) between two 

vectors x and y in C n is given by cos Z(x,y) = [SchOl] . The angle between a vector 

y G C n and a subspace X spanned by I £ C* xp is Z(X,y) = mm x€ x ,|la;||=i ^-(x, y). The 
angle Z(X,Y) between two subspaces spanned by Jf G C" xp and Y G C* xp is defined as the 
largest principal angle between the two subspaces, given by cos /.(X, Y) = a m - m (X H Y) where 
X and Y are orthonormal bases for |_^J and \Y \ j respectively, and cr m ; n denotes the smallest 
singular value |GH06] . The following proposition is a generalization of [AMSV021 Th. 3.1] to 
the complex case. 

Proposition 2.1 Let [X\X±] be a unitary matrix of order n, with X of dimension n x p, 
and let K be an (n — p) x p matrix. Then 

t&nZ(X,X + X ± K) = \\K\\. 

Proof. The matrix Y = (X + X±K)(I + K H K)~ 1 / 2 is an orthonormal matrix with the same 
span as X+X ± K. It follows that cos Z(X, X+X ± K) = a min (X H Y) = a min (I+K H K)~ 1 / 2 = 
(1 + o"max(^0) -1 ^ 2 = (1 + 11-^11 2 ) _1//2 - The conclusion follows from the trigonometric formula 
tan 2 a = (1 — cos 2 a) / cos 2 a. □ 

We now briefly recall some basic facts about invariant subspaces. 

Definition 2.2 (eigenspaces) Let X be a p- dimensional subspace ofC n and let X = [Xi,X2] 
be a unitary nxn matrix such that X\ spans X . Then X H CX may be partitioned in the form 
X H CX = \ ell ) w ^ ere Cu G C pxp . The subspace X is an eigenspace (or invariant sub- 
spacej of C if C21 = 0, i.e., CX C X. By spectrum of X , we mean the set of eigenvalues of 
Cn . We say that X is a nondefective invariant subspace of C if C\\ is nondefective. The in- 
variant subspace X is termed spectral if C\\ and C22 have no eigenvalue in common \GLR86$ . 
The eigenspaces of C H are called left eigenspaces of C . We say that (^L^i?) is a pair of 
spectral left-right eigenspaces of C if and are spectral left and right eigenspaces of C 
with the same spectrum. 

The span of Y G C™ xp is an eigenspace of C if and only if there exists a matrix M such 
that CY = YM. Each spectral eigenspace is isolated, i.e., there exists a ball in Grass(p, n) 
centered on V that does not contain any eigenspace of C other than V. We will also need the 
following result [GV961 §7.6.3], of which we give an informative proof. 

Proposition 2.3 //(3^l,!Vr) is a pair of spectral left-right eigenspaces of C , then there exists 
an invertible matrix S such that the first p columns of S span y^, the first p columns of S~ H 



span y^, and S l CS 



Di 
D 2 



with Di g C pxp . 
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Proof. Let (3^)3^) be a pair of spectral left-right eigenspaces of C. Then there exists X 



unitary such that X CX = q 1 q 1 ^ , the first p columns of X span 3^, and C\\ and C22 
have no eigenvalue in common. Therefore, there exists a matrix L such that C\\L — LC22 = 
— C12 |Ste73] . Let S = X [q j] . Then the first p columns of S span 3^- One easily checks 

Therefore the first p columns 



that S^CS 



Cn 

c 22 



eg 



Moreover C H S~ H = S~ H 

of S~ H span an eigenspace of C H whose spectrum is the same as the one of 3^?- That is, the 
first p columns of S~ H span 3*l- D 
The Rayleigh quotient iteration (RQI) is a classical method for computing a single eigen- 
vector of a Hermitian matrix A. It induces an iteration on the projective space P n_1 that can 
be written as follows. 

Algorithm 2.4 (RQI on projective space) Let A = A H be an n x n matrix. Given Sq 
in the projective space P n_1 , the RQI algorithm produces a sequence of elements o/P n_1 as 
follows. For k = 0, 1, 2, . . ., 

1. Pick y in C n \ {0} such that [y\ = Sk- 

2. Compute the Rayleigh quotient = {y H Ay) / {y H y) . 

3. If A — pkl is singular, then solve for its kernel and stop. Otherwise, solve the system 



(A- Pk I)z = y (1) 



for z. 
4. S k+1 := [z\ 



It is shown in [BS89] that around each (isolated) eigenvector of A, there is a ball in which 
cubic convergence to the eigenvector is uniform. The size of the ball depends on the spacing 
between the eigenvalues. Globally, the RQI converges to an eigenvector for any initial point 
outside a certain set of measure zero described in |BS89| . 

The Grassmann-Rayleigh Quotient Iteration (GRQI) is a generalization of the RQI that 
operates on Grass(p, n), the set of all p-dimensional subspaces of C n |AMSV02] . 

Algorithm 2.5 (GRQI) Let A = A H be an n x n matrix. Given 3^ £ Grass(j>, n), the 
GRQI algorithm produces a sequence of p- dimensional subspaces of C n by iterating from 3\) 
the mapping Grass(j>, n) — > Grass(p, n) : y 1— > 3+ defined as follows. 

1. Pick Y £ CT P such that [Y\ = y. 

2. Solve the Sylvester equation 



AZ - Z(Y H Y)~ l Y H AY = Y (2) 



for Z G C nxp . 
3. Define y + :=[Z\, 



It is shown in [AMSV02] that the subspace 3^+ does not depend on the choice of basis Y for 
3* in the first step. This iteration converges cubically to the p-dimensional eigenspaces of A, 
which are the only fixed points. 

When the matrix A is not normal, the stationary property of the Rayleigh quotient fails. 
Consequently, the convergence rate of the RQI can be at best quadratic. In order to recover 
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cubic convergence, Ostrowski [Ost59b] proposed a two-sided version of the RQI, formulated 
as follows in |Par74j . 

Algorithm 2.6 (Two-Sided RQI) Let C be an n-by-n matrix. Pick initial vectors vq and 
uq satisfying VqUo ^ 0, ||«o|| = ||ito|| = 1- For k = 0, 1, 2, . . 

1. Compute pk = vjfCuk/vjfuk- 

2. If C — pkl is singular solve y H (C — pf.1) = and (C — pkl)x = for y, x ^ and stop, 
otherwise 

3. Solve both vj? +l (C — p^I) = v^u^, (C — pkl)uk+i = u^t^, where and are normalizing 
factors. 

4- U v k+i u k+i = 0; then stop and admit failure. 

The Two-Sided RQI converges with cubic rate to the pairs of left-right eigenvectors of C with 
linear elementary divisor |Par74j . 



3 Two-Sided GRQI 

We propose the following generalization of the Two-Sided RQI, which we call the Two-Sided 
Grassmann-Rayleigh Quotient Iteration (2sGRQI). 

Algorithm 3.1 (2sGRQI) Let C be an n-by-n matrix. Given (3^L >3^io) £ Grass(p, n) x 
Grass(p, n), the 2sGRQI algorithm produces a sequence of pairs of p- dimensional subspaces 
ofC n by iterating from (^Loi^-Ro) the mapping (3^, 3^) ^ (3^+, 3^+) defined as follows. 

1. Pick Yl and Y R in C* xp such that \Yl\ = 3^l and \Yr\ = y^. 

2. Solve the Sylvester equations 

CZ R - Z R {YgY R )- x Y?CY R = Y R (3a) 

V * ' 

Rr 

ZfC - Y»CYr{Y?Yr)- 1 Zf = YE (3b) 



for Z L and Z R in C nx P. 

3. Define y L+ := [Z L \ and y R+ := [Z R \ 



In point 1, one has to choose bases for y^ and y R . There are infinitely many possibilities. 
Indeed, if Y is a basis of y, then {YM : M £ C^, xp } is the (infinite) set of all bases of y. 
Therefore, one has to make sure that 3^l+ an d y R + do not depend on the choice of basis. 
By a straightforward adaptation of the development carried out in [AMSV02J for the GRQI 
algorithm, if (Y L , Y R , Z L , Z R ) solve © then (Y L M, Y R N, Z L M, Z R N) also solve © for all M, 
N in C* xp . Hence, the spans of Zl and Zr only depend on 3 7 l and 3^, and not on the choice 
of the bases Yl and Yr. 

In point 2, the matrix Y^Yr may not be invertible. This corresponds to point 4 in 
the Two-Sided RQI (Algorithm 12. 6p . However, if (3^l ; 3^r) is a pair of spectral left-right 
eigenspaces of C, then Y^Yr is invertible (as a consequence of Proposition I2.3[) . and by 
continuity invertibility holds on a neighborhood of the pair of eigenspaces. 
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In point 2, the (uncoupled) Sylvester equations ([3]) may fail to admit one and only one 
solution. This situation happens if and only if (Yr, Yl) belongs to the set 

S := {(Yl, Y r ) E CT P x C? xp : Rr exists and a(C) n a(R R ) ± 0} 

= (J {(Yl, Yr) G Cr p x Cr p ■ Rr exists and det(R R - XI) = 0}; 

Ae<r(C) 

this follows directly from the characterization of the eigenvalues of Sylvester operators |Ste73l 
Th. 4.4]. Since S is the finite union of algebraic sets, it has measure zero and the interior 
of its closure is empty. This means that if (Yl,Yr) does not yield a unique solution, then 
there exists, arbitrarily close to (Yl,Yr), a pair (Yl,Yr) and a neighborhood of this pair 
on which the solution (Zl,Zr) of ([3|) exists and is unique. In our numerical experiments, 
when such a singularity occurs (i.e., when the solution of the Sylvester equations returned 
by Matlab contains Inf 's or NaN's), we slightly perturb the system. A justification for this 
technique is given in [AMSV02] and the numerical tests performed in Section [8] illustrate that 
the technique works well in practice. 

In point 3, if Zl or Zr is not full rank, then (3^l+, 3^r+) does not belong to Grass(p, n) x 
Grass(p, n). A tall nx p matrix Z is rank deficient if and only if all its p x p minors are zero. 
Therefore, the set 

V := {(Y L ,Y R ) : rank(Zi) < p or rank(Z R ) < p} 

is a subset of a finite union of algebraic sets. So here again, Zl and Zr are full rank for a 
generic choice of Yl, Yr. 

In practice, only a few iterates will be computed. In finite precision arithmetic, the iterates 
no longer improve after a few (typically two or three) iterations because of numerical errors 
(see numerical experiments in Section [8]) . Stopping criteria can rely on the principal angles 
between two successive iterates and on the principal angles between 3^j? and AyR or an d 



4 Practical implementation 

The major computational task in both GRQI (Algorithm 12. 5|) and 2sGRQI (Algorithm 13. 1|) 
is to solve the Sylvester equations. For GRQI, it is recommended to choose an orthonormal 
basis Y (i.e., Y H Y = I p ) that makes Y AY diagonal. This requires solving a p-dimensional 
eigenproblem, which is cheap when p is small. With Y H Y = I and Y AY diagonal, the GRQI 
equation ([2]) decouples into p linear systems for the p columns of Z. We refer to [AMSV02] 
for details. 

The case of 2sGRQI (Algorithm 13. 1 [) is quite different. The matrices Rr and Rl in the 
2sGRQI equations ([3]) are not Hermitian and they may not be diagonalizable. A possible 
approach to solving an equation such as (f3ajl is to reduce it to a certain triangular structure 
by means of unitary transformations and solve the new system of equations using back sub- 
stitution, as described in [GLAM92]. However, we observed in numerical experiments that 
this technique tends to yield rather inaccurate results when the iterates get close to a solution 
(the final error was sometimes around 10 -11 whereas the machine epsilon was approximately 
2.2 • 10~ 16 , to be compared with the results in Table [1]). The reason seems to lie in the fact 
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that the norm of the solution z to the equation (C — pl)z = x becomes very sentitive to p 
when p gets close to an eigenvalue of C. The magic of the classical RQI is that the direction 
of z is well-conditioned, as pointed out by Peters and Wilkinson [PW79] . However, in the 
block case, the magic weakens because, in view of the workings of back substitution, the large 
numerical error in the norm of any one column of Z jeopardizes the accuracy of the other 
columns whose computation depends on that column. 

Consequently, we recommend reducing the small block shifts Rr and Rl in the 2sGRQI 
equations ([3]) to quasi-diagonal form, or to (complex) diagonal form if complex arithmetic 
is available. To fix ideas, assume complex arithmetic and consider the first equation, (|3ap . 
namely 

CZr — ZrRr = Yr. 

Assuming that Rr is nondefective, let Rr = Wr diag(p\, . . . ,p p ) W^ 1 be an eigenvalue de- 
composition of Rr. Multiplying (f3a|) on the right by Wr yields 

CZr - Z R diag(pi, . . . , rhop) = Y R 

where Yr = YrWr and Zr = ZrWr- The advantage of this reformulation of (|3ap is that it 
yields p decoupled equations 

(C - piI)Z R ei = Y r ei 

for each column Zrci of Zr. Back propagation is thus no longer needed. A drawback is 
that this technique does not work when Rr is defective, and numerical errors on Wr may 
become large when Rr is close to being defective. Nevertheless, in our extensive numerical 
experiments on randomly chosen matrices, these difficulties were not noticed (see Section [8]). 

The same kind of discussion applies to the left equation (13bj) . Note that since Rl = 
(Y l h Yr)Rr(Y l h Yr)~ 1 is a similarity transformation of Rr, we have that Wl = (Y l u Yr)Wr 
is a matrix of eigenvector of Rl. Hence, the eigendecomposition of Rl is readily obtained 
from that of Rr. 

5 Local convergence 

The following local convergence analysis can be thought of as a two-sided generalization of the 
proof of cubic convergence of the block-RQI (equivalent to the Grassmann-RQI of |AMSV02j ) 
given in |Smi97j . 

Let (Vl,Vr) be a pair of spectral left-right eigenspaces of C, and let Vl and Vr be cor- 
responding eigenbases. We assume that the eigenspaces are nondefective, that is, the matrix 
(YlVr^^lCVr. is diagonalizable by a similarity transformation. Since (Vi,Vr) is nonde- 
fective, it follows that for all and 3^? sufficiently close to Vl and Vr, the block Rayleigh 
quotients Rr and Rl are diagonalizable by similarity transformations Wr and Wl- Equa- 
tions ([3]) thus can be solved in two steps: (i) diagonalize the small block Rayleigh quotients, 
hence decoupling the equations and reducing them to classical two-sided RQI equations; (ii) 
solve the decoupled two-sided RQI equations, yielding matrices ZlWl and ZrWr that span 
3^l+ and yR+. The key of the convergence analysis is an "oblique" generalization of [SteOlj 
Th. 2], showing that the angles between the right Ritz vectors (the columns of YrWr) and 
the "corresponding" right eigenvectors of A are of the order of the largest principal angle 
between 3^j? and Vr, and likewise for the left Ritz vectors and eigenvectors; see Lemma 15. II 
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Then the result follows quite directly from the cubic convergence of the non-block two-sided 
RQI. 



Lemma 5.1 Let (A, x) be an eigenpair of an nx n matrix C. Let Yr and Yl be orthonormal 
n x p matrices, p < n, such that Y^Yr is invertible. Let wr be an eigenvector of 

B := {Y^YrY^CYr 

associated with the eigenvalue of B that is closest to A. Then 



sin Z(Y R w R ,x) < 



1 + 



2(cos<5) 1 r L as(e) 



sep(w^BwR,w^ ± Bw R± ) - r L ^ 5 (e) 



(1 + tan<5)e 



where e := sin Z(Yr,x) is the angle between the direction of x and the span of Yr, 5 :- 
Z(Yr,Yl) is the largest principal angle between the spans ofYR and Yl, ag(e) :- 



\J1—e 2 — e tan S 

satisfies lim^o ag(e) = 1, 75(e) := (cos<5(Vl — e 2 — etan<5)) _1 (l-|-tan(5)e satisfies lim^o 7<5(e) - 
0, and tl := \\Y£±A h Yl\\ where Yl± € C nx ( n ~ p ) is an orthonormal basis of the orthogonal 
complement of the span of Yl . 

Proof. It is readily checked that the statement is not affected by a unitary change of coor- 
dinates in C n . Therefore, without loss of generality, we work in a unitary coordinate system 



such that Yr 



0(n— p)xp 



. Let Y L ± e C nx ( n " p ) and Y R± £ C nx ( n " p ) be orthonormal bases of 



the orthogonal complements of the spans of Yl and Yr, respectively. Assume without loss of 

X b 

\Y L a 



generality that the eigenvector x has unit norm. Consider the block decompositions x 



and Yl 



Yi 



Lb 



Consider also the decomposition x = Yrxr + Yl±xl±, which yields 



x R := {Y L H Y R )-'Y£x, x l ± := (Y^Y^Y^x. 



-IvH 



H 



-IvH 



Since e = sin Z(Yr,x), we have ||£ a || 2 = 1 — e 2 and ||xb|| = e. We also have (Y^Yr) - ^^ = 
\l T] where T = (Yl ) ^L6- It follows from Proposition 12.11 that ||T|| = tan<5. We also 
obtain 



Yrxr 



[L T] 



x a + Tx b 




Acceptable choices for Yl± and Yr± are Yl± 



l n—p 



(L n _ p + T H T)- l l 2 and Yr± 



px(n—p) 



L n—p 



This yields x L ± = (I n - P + T H T) l / 2 x b and thus ||x£,|| < y/l + tan 2 5e. 
Since sinZ(u, v) < sin Z{u,w) + sin Z(w,v) for all u,v,w € Cq, we have 
Z(Y r wr,x) < Z(Yrwr,Y r xr) + Z(Y R x R ,x). 

Let us first consider the second term in (jl]). Since 
sin Z(Y r xr, x) < \\Y r xr-x\\ < \\Yrxr- 



(4) 



+ 
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it follows that 



smZ(YnXR,x) < \\Txb\\ + \\xb\\ < tan 5 e + e = (1 + tan<5)e. (5) 

Note also for later use that, for all small e such that \/l — e 2 > etan5, we also obtain that 
\\xr\\ > |||a;a|| - \\Txb\\\ > Vl - e 2 - etand. 

We now tackle the first term in Since Yr is orthonormal, it follows that Z(Yrwr, Yrxr) 
Z{wr,xr). Pre-multiplying the equation Cx = \x by (Y l h Yr)~ 1 Y l h yields 

{Y*Y R )- l Y?C{Y R x R + Y L± x L± ) = x R X, 

which can be rewritten as 

[B + E)x R = Xx R , 

where xr := xrHxrH" 1 and 

E := (yf Y R )- l Y? AY^xl^xrW' 1 ^ . 



2\\E\ 



Then, by [JSCM Th. 5.1], 

smZ(w R ,x R ) < tan Z(w R ,x R ) < j- 

sep (w^Bwr, (w R )fB(w R )±) - 2\\E\\ 

if the bound is smaller than 1. The expression of the bound can be simplified using 

\\E\\ = \\(YFY R )- 1 Y?CY L± x L1 .\\\\x R \\- 1 < \\(YtY R )~ l \\\\Y^C H Y L \\\\x L ^Rt l 

< — s r L {l + tan8)e- 



cosi5 y/i - e 2 - etan^' 

where we have used the bound Vl + tan 2 5 < (l+tan<5) that holds for all 5 € [0, f ). Replacing 
all these results in yields the desired bound. □ 

Theorem 5.2 Let (Vl, Vr) be a pair of p- dimensional spectral nondefective left-right eigenspaces 
of an n x n matrix C (Definition \2.2\) . Then there is a neighborhood J\f of (Vl,Vr) in 
Grass(p, n) x Grass(p, n) and a c > such that, for all (VljVr) G M, the subspaces Vl+ and 
Vr+ produced by the 2sGRQI mapping (Alaorithm \3.1\) satisfy 

z(y L+ ,v L ) + z(y R+ ,v R ) < c (z(y L , v L ) + z(y R , v R )f . 

Proof. Since the pair of eigenspaces is assumed to be spectral, it follows that Z(Vl, Vr) < ir/2. 
Therefore, taking the neighborhood J\f sufficiently small, one has Z(Yr, Yl) < 5' < ir/2. More- 
over, since the pair of eigenspaces is assumed to be nondefective, it follows that the eigenbases 
Vr and Vl have full rank. Note that for each column x of Vr, we have Z(Yr, x) < Z(Y R , Vr). 
Lemma I5TT1 implies that for any c\ > 1+tanJ', there exists an e > such that, for all (Vl, Vr) 
with Z(Vl, Vl) + Z(y R , Vr) < e, the angle Z(Yrwr, x) between x and the nearest Ritz vector 
Yrwr satisfies Z(YrWr,x) < c\Z(Yr,x) < c\Z(Yr,Vr). Next, represent the subspaces Vl 
and yR by their Ritz vectors, which decouples (f3a|) into p two-sided RQI equations. By tak- 
ing e sufficiently small, it follows from the cubic convergence of the two-sided RQI that there 
exists C2 > such that, for each column x of Vr, we have Z((zr){,x) < C2(ciZ(Yr,Vr)) 3 for 
at least one column (z R )i of Zr. It follows that Z(Zr,Vr) < c^C2{c\Z{Yr, Vr)) 3 where C3 is 
a constant that depends on the conditioning of the basis Vr. A similar reasoning applies to 
the left subspace. □ 
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6 Comparisons with Newton-based approaches 



It has been long known (see, e.g., Peters and Wilkinson |PW79j ) that the RQI can be viewed 
as a Newton iteration. In fact, the RQI is a Newton method in a certain differential-geometric 
sense [AMS071 . However, the strict interpretation of RQI as a Newton method disappears 
in the block CclS6 ? clS pointed out in [AMSV02J, so much so that the Grassmann-RQI can be 
considered as distinct from the Newton approach. 

Several Newton-based approaches for the general (non-Hermitian) eigenvalue problem 
have been proposed in the litterature. In particular, the well-known Jacobi-Davidson approach 
can be viewed as a Newton method within a sequential subspace algorithm; see, e.g., |LE02t 
IAMS 07] . Here we discuss specifically the Newton method proposed by Chatelin |Cha84j for 
refining eigenspace estimates. The reasoning can be explained as follows. An n x p matrix Y 
spans an eigenspace of C if and only if there exists a p x p matrix M such that 

CY = YM. (6) 

However, any subspace admits infinitely many bases, and the solutions Y of ([6|) are thus not 
isolated. A way to remove the freedom in the choice of basis is to impose on Y a normalization 
condition W H Y = I where W is a given full-rank n x p matrix. Then ([6]) becomes 

F(Y) := CY - Y(W H CY) = (7) 

where the unknown Y is normalized by W H Y = I. The Newton iteration for solving ([7J) is 
given by 

(7 - YW H )CA - A(W H CY) = —F(Y), W H A = (8) 
Y + := Y + A. (9) 

If the basis Y is chosen orthonormal and W := Y, then ([H]) becomes 

nCTIA - A{Y H CY) = -nCY, Y H A = (10) 

where H := I — YY H . The resulting algorithm admits an interpretation as a Newton method 
on the Grassmann manifold }A~MS07j . The rate of convergence is quadratic in general (cubic 
when C is Hermitian). 

The constraint Y H A = can be addressed by setting A = Y±K, where Y± is an or- 
thonormal matrix with Y H Y± = and K is an (n — p) x p matrix; see, e.g., |Dem87j . Then 
Y H A = is trivially satisfied and equation (I10p becomes 

{YfCY ± )K - K(Y H CY) = -YfCY, (11) 

i.e., a Sylvester equation without constraints on the unknown K. As pointed out in |AMS V02] . 
solving (llip takes 0(ra 3 ) operations even when C is condensed (e.g. tridiagonal) because 
Y^CY± is a large dense (n — p) x (n — p) matrix. However, Lundstrom and Elden proposed 
an algorithm |LE021 alg. 2] for solving ()10p that does not require the computation of Y±CY±. 
It takes 0(np 2 ) operations to solve (|10p when C is block diagonal of sufficiently moderate 
block size and 0(n 2 p) when C is Hessenberg. The complexity of the 2sGRQI method (Algo- 
rithm 13. ip is of the same order. 
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A theoretical comparison between algorithms based on inverse iteration and on Newton 
does not reveal that one approach has a clear edge over the other. Among the advantages 
of the 2sGRQI method (Algorithm 13. 1 H over Chatelin's method, one can mention that the 
convergence of 2sGRQI is cubic instead of quadratic, and that a pair of left-right eigenspaces is 
computed instead of just a right-eigenspace. On the other hand, Chatelin's method admits a 
convergence analysis when the target eigenspace is defective |Dem87[ AMS04J . and it requires 
solving only one Sylvester equation instead of two in 2sGRQI. However, we show in Section [7] 
that one Sylvester equation suffices for 2sGRQI on some important structured eigenproblems. 



7 Structured eigenproblems 

In this section, we show that the 2sGRQI induces particular one-sided formulations for some 
structured eigenproblems. 



7.1 .E-Hermitian eigenproblem 

Let C be an n x n matrix. If there exists an invertible matrix E such that 

EC = C H E, (12) 

then we say that C is E-Hermitian. If C is .E-Hermitian, then its left and right eigenspaces 
are related by the action of E. Indeed, let S be a (complex) matrix of principal vectors of C, 
i.e., 

CS = SD 

where D is a (complex) Jordan matrix; then, from (I12p . one obtains C H (ES) = (ES)D. 

The case where E is Hermitian or skew-Hermitian, i.e., E H = ±E, is of particular interest 
because, as we show in the next proposition, the relation = Ey^ is invariant under the 
2sGRQI (Algorithm 13 . X [> . Therefore, if = EyR, it is not necessary to solve both (f3a|) 
and (|3bj) : just solve (f3a|) to get 3^+, and obtain 3^+ as yL+ '■= EyR + . Moreover, since the 
pairs of left-right eigenspaces of C also satisfy Vl = EVr, Theorem 15.21 also applies. 



Proposition 7.1 Let E be invertible with E H = ±E and let C be E-Hermitian, i.e., EC 
C H E. If Yl = EY, Y R = Y, and Z satisfies 



CZ-Z (Y H EY)~ 1 (Y H ECY) = Y, 



(13) 



then Zl = EZ and Zr = Z satisfy the 2sGRQI equations Q. Hence, if y^ = Ey^, then 
3^l+ = Ey^ + . Moreover, the subspace iteration [Y\ \— > \ Z\ defined by (|13|) converges locally 
cubically to the spectral nondefective right- eigenspaces of C. 



Proof. It is easy to check that replacing Yr := Y, Zr := Z, Yl := EYr, Zl := EZr in (|3ap 
and (pb|) yields (fTHj) in both cases. In order to prove cubic convergence, it is sufficient to notice 
that the pairs (Vl,Vr) of eigenspaces satisfy Vl = EVr, as was shown above. Therefore, if 
y is close to Vr, then the pair (3^,3^) := (Ey,y) is close to {Vl,Vr) and local cubic 
convergence to Vr_ follows from Theorem 15.21 □ 

The discussion in Section H] on solving Sylvester equations applies likewise to (1131) . 
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Generalized Hermitian eigenproblem 



Using Proposition l7.ll we show that the 2sGRQI yields a Grassmannian RQI for the Hermitian 
generalized eigenproblem AV C BV which does not involve an explicit computation of B~ 1 A. 
Let A and B be two Hermitian n-by-n matrices with B invertible. Consider the problem of 
finding a p-dimensional subspace V such that AV C BV. Let V G <C nxp be a basis for V, then 
AV C BV if and only if there is a matrix M such that AV = BVM. Equivalently, V spans 
a right-eigenspace of B _1 A, i.e., 

B~ X AV = VM. 

The problem is thus to find a right-eigenspace of C := B _1 A. The conditions in Proposi- 
tion [T7TJ are satisfied with £7 := B. The modified GRQI equation (| 13j) becomes 



- sz (y^^y)" 1 (y H ^y) = #y (14) 



and the subspace iteration [Y\ i— ► |_ZJ converges locally cubically to the spectral nondefective 
eigenspaces of B _1 A. In particular, S" 1 ^! is nondefective when A or i? is positive definite. 



Skew-Hamiltonian eigenproblem 



Let T be a skew-Hamiltonian matrix, i.e., (TJ) H = —T J, where J = ( q ) , see e.g. [BBMX02] , 
Equivalently, JT = T H J, i.e., T is J-Hermitian. Conditions in Proposition 17.11 are satisfied 
with C := T and E := J. The modified GRQI equation (jl3l) becomes 



rz - z (y H jy)" 1 (y H jty) = y 



(15) 



and the subspace iteration \Y\ ^ \_Z\ converges locally cubically to the spectral nondefective 
right-eigenspaces of T. 



7.2 E'-skew-Hermitian eigenproblem 

Let E be an invertible n x n matrix and let C be an E- skew- Hermitian nx n matrix, namely 

EC = -C H E. (16) 

We saw in the previous section that the corresponding left and right eigenspaces of E- 
Hermitian matrices are related by a multiplication by E. The case of i^-skew-Hermitian 
matrices is slightly different. 



Proposition 7.2 Let C be an E- skew- Hermitian matrix. Then the spectrum of C is sym- 
metric with respect to the imaginary axis. In other words, if A is an eigenvalue of C , then 
so is —A. Moreover, ifVi and Vr are left and right eigenspaces of C whose spectra are the 
symmetric image one of the other with respect to the imaginary axis, then Vl = EVr. 

Proof. Letting S be an invertible matrix of principal vectors of C, i.e., 

CS = SD (17) 
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where D is a Jordan matrix, (|16|) yields 

C H ES = ES(-D). (18) 

Hence, the matrix —D is a Jordan matrix of C . Therefore, if A is an eigenvalue of C , then 
— A is an eigenvalue of C , and thus —A is an eigenvalue of C. Moreover, equations (|17p 
and (fT8|) show that if V is a right-eigenspace of C with eigenvalues Ai x , . . . , Aj p , then EV is a 
left-eigenspace of C with eigenvalues — Aj 1} . . . , — Aj . □ 

Consequently, letting V be a spectral right-eigenspace of C, we have that (EV, V) forms a 
pair of spectral left-right eigenspaces of C if and only if the spectrum of V is symmetric with 
respect to the imaginary axis. We call such an invariant subspace V a full eigenspace of the 
E'-skew-Hermitian matrix C. 

If E is Hermitian or skew-Hermitian, then the relation y^ = Ey R is invariant by the 
2sGRQI (Algorithm 13. IB . as we show in the forthcoming proposition. Therefore, if y^ = Ey R , 
it is sufficient to solve (j3a|) only, and then compute 3^l+ := Ey R+ . Moreover, the 2sGRQI 
iteration restricted to the pairs (3^,3^) = (Ey,y) converges locally cubically to the full 
nondefective eigenspaces of C. 



Proposition 7.3 Let E be invertible with E H = ±E and let C be E -skew-Hermitian, i.e., 
EC = -C H E. If Y L = EY and Y R = Y, then Z L = -EZ and Z R = Z satisfy the 2sGRQI 
equations ([3]) with 

(19) 



CZ - Z (Y H EY)~ 1 (Y H ECY) = Y. 



Therefore, ify L = Ey R , then y L+ = Ey R+ . 

Moreover, let V be a full nondefective right-eigenspace of C (which means that the eigenvalues 
of C|v have the same multiplicity as in C , the spectrum of C|y is symmetric with respect to 
the imaginary axis, and C|v is nondefective). Then the subspace iteration \Y\ ^ [Z\ defined 
by (|13p converges locally cubically to V. 



Note that this proposition differs from Proposition 17.11 in two points: Zl = —EZ and the 
specification that V must be full. 

Proof. It is easy to check that replacing Y R := Y, Z R := Z, Yl := EY R , Zl := —EZ R in (f3"a|) 
and (|3bp yields (TL9j) in both cases. In order to prove cubic convergence, it is sufficient to 
notice that the pairs (Vl,V r ) of full nondefective left-right eigenspaces satisfy Vl = EV R , as 
was shown above. Therefore, if y is close to V R , then the pair (3^,3^) := (Ey,y) is close 
to (Vl, V r ) and local cubic convergence to V follows from Theorem 15.21 □ 



Skew-Hermitian eigenproblem 



Let fl be skew-Hermitian. Then we have EC 
modified GRQI equation (I19p becomes 



C H E with C := tt and E := I. The 



nZ - Z (Y H Y)~ 1 (Y H nY) = Y. 



(20) 



This is simply the classical GRQI equation - This is not surprising as skew-Hermitian 
matrices are normal matrices. 
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Hamiltonian eigenproblem 



Let H be Hamiltonian, i.e., (HJ) H = HJ, where J = (_?jg). 
-H H J. Thus we have EC = -C H E with C := H and E := 
equation (fTUj) reads 



HZ - Z (Y H JYy 1 (Y H JHY) = Y. 



This is equivalent to J-ff = 
J, and the modified GRQI 

(21) 



Proposition [73] implies that the subspace iteration with iteration mapping [^J •— > \_Z\ defined 
by (I2ip converges locally cubically to the full nondefective right-eigenspaces of H. 



7.3 The generalized eigenvalue problem 

We briefly discuss the application of the 2sGRQI concept to the generalized eigenvalue prob- 
lem. Let A, B S C nxn . The generalized eigenvalue problem consists in finding the nontrivial 
solutions of the equation Ax = XBx. Corresponding to the notion of invariant subspace 
for a single matrix, we have the notion of a deflating subspace, see e.g. |Ste73[ IGV96j . The 
p-dimensional subspace X is deflating for the pencil A — XB if there exists a p-dimensional 
subspace y such that 

AX,BX C y. (22) 

Here we suppose that the pencil A — XB is nondegenerate, i.e., det(A — XB) is not trivially 
zero. Then there exists a and f3 such that B := aB — (3A is invertible. Now take 7, 5 such 
that a8 - 7/3 7^ and let A := jB - 5 A. Then ([22]) is equivalent to 

B^AX C X 
BX = y, 

i.e., X is an invariant subspace of B~ 1 A. Replacing this expression for C in ([3]), one obtains 
after some manipulations 

AZ R Y^BY R - BZ R Y^AY R = BY R (23a) 
A H Z L Y^B H Y L - B H Z L Y§A H Y L = B H Y L (23b) 

where Yl '■= B~ h Yl and Zl := B~ H Z^. It yields an iteration for which Yr and Yl locally 
cubically converge to pairs of left-right deflating subspaces of the pencil A — XB. Note that 
if B is invertible then we can choose B := B and A := A. 



8 Numerical experiments 

We report on numerical experiments that illustrate the potential of the 2sGRQI method 
(Algorithm 13. 1|) as a numerical algorithm. The 2sGRQI method has been implemented in 
Matlab as described below. 

Algorithm 8.1 (implementation of 2sGRQI) Let C be annxn matrix. Given two nxp 
matrices Yl and Yr satisfying Y^Yl q = I = Y^Yr , the algorithm produces a sequence of 
matrices {YL k > ^R k ) as follows. For k = 0, 1, 2, ... , 

1. Compute the p x p block Rayleigh quotient Rr := {Y^YR k )~ 1 Y[^CYR k . Compute an 
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eigendecomposition Rr = Wijdiag(pi, . . . , p p )W R using the Matlab eig function. Obtain the 
eigendecomposition Rl = Wj^diag(pi, . . . ,p p )W£ H by computing Wl '■= (Y^Yr^Wr. 

2. Solve the decoupled equations ((3|), that is, (C — Pil){zjt)i = Y^WrCi and (C H — Pil){zi)i = 
Yi k W£ H Ci, i = 1, . . . ,p, using the Matlab '\" operator. If the solutions have any nonfinite 
element, then solve instead (C — pil + el)(zjt)i = YR k WRei and (C — pil + el){zi)i = 
Yi k W£ H ei, i = 1, . . . ,p, with e small (we took e = 10 3 u||C||,f where u is the floating point 
relative accuracy and \\C\\f is the Frobenius norm ofC). 

3. Orthonormalize Zr := \(zr)i ••• (<zr)i] to obtain Yji k+X , and likewise for Zr to obtain 
Y^ k . In Matlab, orthonormalizations are performed using the "economy size" QR decompo- 
sition, [YL, ignore] = qr(ZL,0) and [YR, ignore] = qr(ZR,0). 

Note that if C, Yl and Yr are real, then the columns of Zl and Zr appear in complex 
conjugate pairs and unnecessary work can thus be avoided in the computation of Zl and Zr. 

It is well known [BS89] that the basins of attraction of RQI (Algorithm I2.4p may collapse 
around attractors when the eigenvalues of A are not well separated. This property also 
holds for GRQI |ASVM04| and obviously extends to 2sGRQI (Algorithm EU) . Moreover, in 
2sGRQI the matrix C is not necessarily Hermitian; its eigenspaces can thus be arbitrarily 
close to each other. In a first set of experiments, in order to ensure a reasonably large basin 
of attraction around the left-right eigenspaces, we ruled out clustered eigenvalues and ill- 
separated eigenvectors by choosing C as follows: C = SDS' 1 , where D is a diagonal matrix 
whose diagonal elements are random permutations of 1, . . . ,n and S = I + »^» E, where 
the elements of E are observations of independent random variables with standard normal 
distribution and a is chosen from the uniform distribution on the interval (0, 0.1). The initial 
matrices Yl q and Yr are randomly chosen such that dist([5 / R J, [S(:,l : p)\) < 0.1 and 
dist(|_iz, J , [S~ H (:, 1 : p)\) < 0.1, where "dist" is the largest principal angle. 

Algorithm 18.11 was run 10 6 times with n = 20, p = 5. The matrices C, Yl , and Yr 
were randomly chosen in each experiment as explained above. Experiments were run using 
Matlab 7.2 with floating point relative accuracy approximately equal to 2 • 10~ 16 . Results are 
summarized in Table [TJ where the error e is defined as the largest principal angle between 
[Yr\ and [S(:, 1 : p)\ plus the largest principal angle between [Yl\ and [S~ H (:, 1 : p)\. These 
results show that convergence to the target eigenspace occurred in each of the 10 6 runs. The 
evolution of the error is compatible with cubic order of convergence. 



Iterate number 


mean(logl0(e)) 


max(logl0(e)) 





-1.4338 


-1.0000 


1 


-4.6531 


-2.6338 


2 


-13.9359 


-8.3053 


3 


-16.5507 


-15.1861 


4 


-16.5524 


-15.1651 


5 


-16.5509 


-15.1691 



Table 1: Numerical experiments for Algorithm 18.11 See details in the text. 

The behavior of the 2sGRQI algorithm in case of ill-separated eigenvectors/values would 
deserve investigation. The Hermitian case is studied in [ASVM04] where improvements of 
GRQI and the Riemannian Newton algorithm are proposed. 
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In another set of experiments, real Hamiltonian matrices C were selected randomly as 



C 



F 

H + H H 



G + G> 
—F H 



where F, G and H are matrices of dimension | x | whose elements are independent obser- 
vations of the standard normally distributed random variable. A new matrix C was selected 
for each experiment. For testing purposes, an eigenvalue decomposition C = SDS" 1 was 
computed using the Matlab eig function, and the full left and right real eigenspaces corre- 
sponding to the eigenvalues with largest real part in magnitude were chosen as the target 
left and right eigenspaces. (The notion of full eigenspace is defined in Section 17.21 The 
real eigenspace associated to a pair (A, A) of complex conjugate eigenvalues with eigenvectors 
v r + ivi and v r — ivi is the span of v r and Vi.) The eigenvalue decomposition was ordered 
in such a way that [S~ H (:,l : p\) is the target left-eigspace and [S" (: , 1 : p)\ is the target 
right-eigenspace. Note that we have p = 2 when the target eigenvalues are real (A and —A), 
or p = 4 when the target eigenvalues have a nonzero imaginary part (A, A, —A, and —A). The 
initial matrix Y# was randomly chosen such that dist ( [y/? J , [S(-, 1 : p)\) < 0.1, where "dist" 
is the largest principal angle, and Yl was chosen as JYr in accordance with the material 
of Section 17.21 Convergence to the target left and right eigenspaces was declared when the 
error e as defined above was smaller than 10 -12 at the 10th iterate. Algorithm 18.11 was run 
10 6 times with n = 20 and p = 4 with the matrices C, Y^ and Yr randomly chosen in each 
experiment as described above. Note that, in accordance with the material in Section 17.21 
only Zr was computed at each iteration; Zl was chosen as JZr. We observed that conver- 
gence to the target eigenspaces was declared for 99.95% of the 10 6 experiments. Next, the 
experiment was run 10 6 times with the distance bound on the initial condition set to 0.001 
instead of 0.1. Convergence to the target eigenspaces was declared for all but seven of the 10 6 
randomly generated experiments. This confirms the potential of Algorithm 18.11 for refining 
initial estimates of full eigenspaces of Hamiltonian matrices. 

9 Conclusion 

We have shown that Ostrowski's two-sided iteration generalizes to an iteration on Grass(p, n) x 
Grass(p, n) that converges locally cubically to the pairs of spectral nondefective left-right 
eigenspaces of arbitrary square matrices. The iteration is competitive with Chatelin's New- 
ton method and it yields one-sided formulations adapted to some structured eigenproblems, 
including the Hamiltonian and generalized Hermitian eigenproblems. 
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